## load necessary libraries
##
library(Cairo)

## initialize device
##
#Cairo(file="rademacher.png", type="png", width = 7, height = 7, pointsize = 12, record = getOption("graphics.recorecord"),
# rescale = c("R", "fit", "fixed"), xpinch, ypinch, bg = "transparent", canvas = "white",
# gamma = getOption("gamma"), xpos = NA, ypos = NA, buffered = getOption("windowsBuffered"),
# restoreConsole = TRUE)
CairoPNG(filename = "Rademacher%03d.png", width = 760, height = 350, pointsize = 12, bg = "white")
par(mfrow=c(1,2))

## Rademacher func definition
##
rademacher <- function(x, ofst, amp, per){
 amp*sign(sin( 2*pi*(x-ofst) / per )) 
}

## Some global variables
##
sampling.frequency <- 1024 
sampling.delta <- 1/sampling.frequency 
rademacher.offset <- 0.0
rademacher.period <- 0.2
rademacher.amplitude <- 0.8

## Plot rademancher func
##
sample.t<-seq(-0.5,1.5,by=sampling.delta)
sample.v<-rademacher(sample.t, ofst=rademacher.offset, amp=rademacher.amplitude, per=rademacher.period)
plot(0,0,col="#00000000",ylim=c(-1,1), xlim=c(-0.5,1.5),xlab="X",ylab="Y",main="Rademacher function, 5Hz")
lines(sample.t, sample.v, lwd=2, col=heat.colors(32)[15])
abline(h=0, lwd=1, col=heat.colors(32)[1])
abline(v=0, lwd=1, col=heat.colors(32)[1])

## FFT stuff goes in here
##
sample.fft <- fft(sample.v)
sample.f.Nyquist <- 1 / 2 / sampling.delta
sample.f <- sample.f.Nyquist * c(seq(length(sample.t)/2), -rev(seq(length(sample.t)/2))) / (length(sample.t)/2)
## plot it
plot(sample.f, Mod(sample.fft)[2:length(Mod(sample.fft))]/length(Mod(sample.fft)-1), xlim=c(-10,10), type='o', lwd=2, 
 col=heat.colors(32)[15], xlab="Frequency, Hertz", ylab="Power", main="Simple spectral analysis")

dev.off()




##
##
## The end, code below is for development :)
##
##
##
##





## Convert into time series (skip that) and get default spectrum function to work
##
#spectrum(y, method="ar", main="R embedded spectrum() func", sub="method=\"ar\"")
#spectrum(y, method="pgram", main="R embedded spectrum() func", sub="method=\"pgram\"")

## FFT time in here
##
# 1.0 prepearing values vector for fft to be power of 2 length with zero padding 
pow <- 0
while( length(sample.v) > 2^pow ) pow <- pow+1 
sample.nfft <- 2^pow
fft.values <- rep(0.0, sample.nfft)
fft.values[1:length(sample.v)] <- sample.v
## Getting fft for our sample
sample.fft <- fft(fft.values)
## FFT is symmetric - get rid of second half () and consider only DC part by taking Modulus which == ABS
sample.fft <- Mod(sample.fft[1:ceiling((sample.nfft+1)/2)])
sample.fft.scaled <- sample.fft/length(sample.t)
## Take the square of the magnitude of fft of x.
sample.fft.scaled <- sample.fft.scaled^2
## Since we dropped half the FFT, we multiply mx by 2 to keep the same energy.
## The DC component and Nyquist component, if it exists, are unique and should not
## be mulitplied by 2.
if ( length(sample.fft.scaled) %% 2 == 0){
 sample.fft.scaled[2:(length(sample.fft.scaled)-1)] = sample.fft.scaled[2:(length(sample.fft.scaled)-1)]*2;
}else{
 # odd nfft excludes Nyquist point
 sample.fft.scaled[2:length(sample.fft.scaled)] = sample.fft.scaled[2:length(sample.fft.scaled)]*2;
}
sample.freq.scale = rep(0:(length(sample.fft.scaled)-1),by=1)*(sampling.frequency/sample.nfft)
plot(sample.fft.scaled, col="red", type="o", xlim=c(0,10), xlab="Fourier Frequency")

################
################

## Extract DC component from transform 
dc <- Mod(sample.fft)/N
periodogram <- round( Mod(sample.fft)^2/N, 3)
periodogram <- periodogram[-1]
periodogram <- periodogram[1:(N/2)]
maxfreq <- s$freq[ which(max(s$spec) == s$spec) ]
print(1/maxfreq)
plot(periodogram, log='y', main='FFT of ECG vs frequency', xlab='Frequency [Hz]')
#
#s <- spectrum(y, taper=0, detrend=FALSE, col="red", main="Spectral Density") 
#plot(log(s$spec) ~ s$freq, col="red", type="o", xlab="Fourier Frequency", xlim=c(0.0, 0.5), ylab="Log(Periodogram)", main="Periodogram from 'spectrum'")


sample.total.num <- length(sample.v) # total number of the samples
sample.total.time <- sample.t[length(sample.t)]-sample.t[1]
delta.time.between.samples <- samples.total.time / samples.total.num
sampling.frequency <- 1 / delta.time.between.samples
## The Nyquist frequency (folding frequency, or the cut-off frequency) which is is half the sampling frequency
sampling.nyquist <- 1/2/delta.time.between.samples
## Calculating frequencies points
sample.frequencies <- sampling.nyquist * c(seq(length(y)/2), -rev(seq(length(y)/2))) / (length(y)/2)




set.fft<-fft(y)
## Plot fft():
##
#plot(set.fft, type='l',lwd=1, col=heat.colors(32)[15])
## Plot Mod(fft()):
##
plot(Mod(set.fft), type='l', log='y', main='FFT of vs index',lwd=1, col=heat.colors(32)[15]) 

## Now, sample period is 1.0
delta <- 1.0
sampling.freq <- 1.0/sampling.delta